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ABSTRACT 

The issue of fragmentation in self-gravitating gaseous accretion discs has implications 
both for the formation of stars in discs in the nuclei of active galaxies, and for the 
formation of gaseous planets or brown dwarfs in circumstellar discs. It is now well 
established that fragmentation occurs if the disc is cooled on a timescale smaller than 
the local dynamical timescale, while for longer cooling times the disc reaches a quasi- 
steady state in thermal equilibrium, with the cooling rate balanced by the heating 
due to gravitational stresses. We investigate here how the fragmentation boundary 
depends on the assumed equation of state. We find that the cooling time required 
for fragmentation increases as the specific heat ratio 7 decreases, exceeding the local 
dynamical timescale for 7 = 7/5. This result can be easily interpreted as a consequence 
of there being a maximum stress (in units of the local disc pressure) that can be 
sustained by a self-gravitating disc in quasi-equilibrium. Fragmentation occurs if the 
cooling time is such that the stress required to reach thermal equilibrium exceeds this 
value, independent of 7. This result suggest that a quasi-steady, self-gravitating disc 
can never produce a stress that results in the viscous a parameter exceeding ~ 0.06. 
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1 INTRODUCTION 

It is becoming clearer that self-gravity may play an 
important role in the dynamical evolution of accretion 
discs. Active galactic nuclei (AGN) often show rota- 
tion cu rves that depart s i gnific a ntly from a Keplerian 
profile (iGreenhill fc Gwinnl Il997fc lLodato fc Bertinl l20od 
iKondratko et al-lbOOoTT while protostellar disc masses may 
be a significant fraction of the c entral object mass during the 
early stages of star formation <Lin fc Pringlelll987l : I Larson 
ll984Tl . For example, a recent radio observation of a Class 
object (a very young prot ostellar source, with a ge < 10 5 
yrs and mass M* m 0.8Mm. lRodriguez et al.ll2005l) shows a 
disc with a estimated mass Mdi ac = 0.3 — OAMq and a very 
suggestive two armed spiral-like stucture. In addition, there 
are also clues that massive accretion discs can be present 
aroun d massive protostars feeltran et al.ll2u04t IChini et alJ 
12004) . 

An important effect of disc self-gravity is that it pro- 
vides an efficient mechanism for transporting angular mo- 
mentum out wards, allowing mass to accrete onto the cen- 
tral object llLin fc Pringld Il987l : lLaughlin fc BodenheimeU 
1994). In protostellar discs, in which the ionization level 



is exp ected to be low, thus inhib i ting MHD-driv en turbu- 
lence jMatsumoto fc Taiimalll995l : lGammiellT99r3) . this may 
be the dominant mechanism, at least in certain regions of 
the disc. In general, however, the transport associated with 
self-gravitating disturbances may not be well described as 
a simple diffusion mechanism. iBalbus fc Papaloizoul jl999l) 
have shown that, for self-gravitating discs, the energy flux 
contains non-local terms that they a ssociate with wave en- 
ergy transport. On the other hand. lLodato fc Ricd (I2004I . 
2005) have shown that, for a self-gravitating accretion disc 
in thermal equilibrium, the dissipation arising from gravita- 
tional stresses agrees reasonably well with the expectations 
based on the stand ard viscous theo ry <|S hakura fc Sunvaevl 
Il973t lPringlelll98ll) (see, however. iMeiia et al.ll2005l for a 
different result). 

I n the simulations performed bv lLodato fc Rice! (I2004I . 
2005) thermal equilibrium is achieved by allowing the disc 
to heat up through PdV work and shock dissipation, and 
cooling it at a prescribed rate. They use a simple cooling 
term with a cooling time given by t coo i — 1 , where Q is 
the angular frequ ency and = 7.5. In thermal equilibrium 
it can be shown (|PringleMl98lT: lGammiell200lD that, for a 
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viscous disc, the parameter a iShakura fc Sunvaevlll97sl) . 
that characterises angular momentum transport, and the 
cooling time, t coo i, are related through 

4 1 

97(7- l)tcoolfi' 1 ' 

where 7 is the ratio of the specific heats. 

The ultimate evolution of a self-gravitating accretion 
disc depends strongly on the rate at which the disc heats 
up, through the growth of the instability, and on the rate 
at which it cools. I t has been sho wn, using a local, two- 
dimensional model jGammiellioOlh . that if i coo i < 3fl , 
the disc will fragment into bound objects rather than evolve 
int o a quasi-st e ady s tate, a result that was largely confirmed 
bv lRice et al.l ((2003) using global, three-dimensional mod- 
els. This process has been suggested as a mechan ism for 
forming both ga s eous planets in pr otostellar discs IlKuipeil 
1 195 it lBosdll998t iMaver et all 12004 . and for f orming stars 
in active galactic discs IIGoodman fc Tanl l2004), 

A major obstacle to the formation of objects via 
disc fragmentation is the requirement that th e cooling 
time be smaller than the local dynamical time iGammid 
l200ll:lRice et al]l2003ft. Th is may be possible in AGN disc 
jjohnson fc Gammic 2003J), but appears unlikely in proto- 
stellar discs (lRafik^vll2005r) . The cooling time requirement 
was, however, determined us ing equations o f state with spe- 
cific heat ratios of 7 = 2 jGammiell200ll), and 7 = 5/3 
jRice et alJ 12003ft . It has been suggested llLodato fc Ricd 
l2005ft" ~t"hat since the stress required to balance the cooling 
rate (as measured by the viscous a) depends on the specific 
heat ratio, 7 (see equation 0, the cooling time required for 
fragmentation may also depend on 7. In particular, if there 
is a maximum stress a max (in units of the local disc pres- 
sure) sustainable by a self-gravitating disc, then equation 
Q states that fragmentation should occur if: 



^cool^ ^ 
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(2) 



In this paper we investigate, using global, three- 
dimensional simulations, how the fragmentation boundary 
(as measured by t coo i) varies for different values of 7. We 
also consider var ious disc masses to study the suggestion by 
iRice et all J2003h that the fragmentation boundary may de- 
pend on the ratio of the disc mass to the mass of the central 
object. In §2, we describe the range of cooling times, disc 
masses, and specific heat ratios, 7, that have been consid- 
ered, and determine, for a given disc mass, the cooling time 
required for fragmentation. In §3 we discuss these results in 
light of the relationship between the stresses in the disc (as 
measured by the viscous a) and the imposed cooling time. 
We show that fragmentation is indeed easier in discs with 
smaller specific heat ratios. We therefore conclude that there 
is a maximum stress sustainable by a self-gravitating disc 
and we quantify this maximum stress to be a max ~ 0.06. In 
§4 we discuss our results and draw our conclusions. 



2 SIMULATION RESULTS 

Th e simulatio n s perf ormed here are very s i milar t o thos e 
of IRice et alJ l|2003|) and lLodato fc Rice! (|2004l . l200it . 

The three-dimensional, gaseous discs are modelled using 



Smoothed Particle Hydrodynamics (SPH) jMonaghanll992l ; 
lBenzlll99Cft . a Lagrangian hydrodynamics code. The disc is 
represented by 250000 SPH particles, while the central star 
is a point mass onto which gas particles may accrete if they 
approach to within the accretion radius (here taken to be at 
a radius of R acc = 0.25). In code units, the disc extends from 
Rin = 0.25 to i?out = 25, and the central object has a mass 
of M* = 1. We consider disc masses of M^isc = 0.1, 0.25 and 
0.5, with initial surface density profiles of E oc r _1 , and with 
a temperature that has a radial profile of T oc r~ 0,5 . With 
these initial surface density and temperature profiles, the 
Toomre stability parameter, Q = c 3 Q/nG~S, is not initially 
constant, but decreases with increasing radius. The temper- 
ature is therefore normalised such that at the beginning of 
the simulation the disc is stable, with a minimum Q = 2 at 
R = 25. 

Since we are interested in how cooling influences the 
disc evolution, the disc gas is allowed to heat up due to 
both PdV work and viscous dissipation, with t he viscosity 
given by the standard SPH artificial viscosity iMonaehanl 
1992) with qsph = 0.1, and /Ssph = 0.2. We use an adiabatic 
equation of state and consider specific heat ratios of 7 = 5/3, 
and 7 = 7/5. A particle with internal energy per unit mass 
Ui is then cooled using 



dui 
~dT 



Ui 
tcooi 



(3) 



where t coo i = /3f2 _1 . 

An important numerical issue to be considered in this 
context is the role of artificial SPH viscosity. The growth 
and the saturation of gravitational instabilties depends on 
the balance between external cooling and internal heating 
provided by the instability itself. Therefore, we have to be 
sure that dissipation is dominated by gravitational instabil- 
ities rather than by artificial viscosity (which would provide 
an extra, undesired, s tabilization term in the energy bal- 
ance) . It can be shown iArtvmowicz fc Lubowll993lMurravl 
1996) that artificial SPH viscosity scales as vsph oc qsph (ft), 
where (h) is the average SPH smoothing length. We there- 
fore should achieve a high resolution (in order to keep the 
smoothing length small) and adopt a sufficiently low value 
of osph, while preserving the ability of the code to prop- 
erly resolve the shocks that arise in th e simulation. We have 
already shown (lLodato fc Ricdl2004l Appendix) that with 
the setup described above (N = 250000, qsph = 0.1 and 
/3 SPH = 0.2), we are indeed able to properly resolve the 
shocks and to have an artificial dissipation smaller by more 
than one order of magnitude with respect to gravitationally 
induced dissipation. We are therefore confident that artificial 
viscosity is not going to affect significantly our conlcusions. 

For each disc mass, and for each 7, we have performed 
a large number of simulations with different values of /3. We 
initiall y start with a value that shoul d result in fragmen- 
tation iGammidl200lt IRice et alJl2003Tl . We stop the sim- 
ulation once at least one clump/fragment has formed that 
has a density 2 — 3 orders of magnitude greater than the 
surrounding gas. The densest clump is then tested to check 
if it is bound. Firstly, we determine the approximate size of 
the clump, by finding the distance from the center of the 
clump at which the density has returned to a value com- 
parable with that of the surrounding disk. All the particles 
within this spherical volume are then assumed to be part of 
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Figure 1. Surface density structure of discs with masses M^ ac = 
0.1 and with cooling times of t coo \Q = 3 (top left), t coo if2 = 5 
(top right), t coo if2 = 6 (bottom left), and t coo iQ = 7 (bottom 
right). The logarithmic colour scale in each figure is from 10 g 
cm -2 to 2 X 10 4 g cm -2 . The linear scale is from -25 to 25 for 
both axes. 



this clump. In every case in which fragmentation occurred, 
the densest clump consisted of at least 100 particles, and in 
some cases as many as 500 particles. This mo re than satis- 
fies the Jeans criterion feate fc Burkert)ll997l) . and we are 
therefore confident that the fragmentation in these simula- 
tions is not artificial. Once the clump size has been deter- 
mined, we then calculate the total thermal energy and the 
gravitational potential energy. If the net energy is negative 
the clump is bound, the simulation is stopped, and a new 
simulation is started with the same initial conditions, but 
with /3ncw = /Sold + 1. This new simulation is then run until 
it either also frag ments or until the disk has evolv ed into a 
quasi-steady state <lGammiell200 it iRice et ai1l2003l) without 
any signs of fragmentation. We also run the non-fragmenting 
simulations approximately an outer rotation period longer 
than the equivalent simulation that did undergo fragmenta- 
tion. 

Figure Q illustrates the procedure discussed above. The 
disc mass in all four figures is Mdisc = 0.1, the specific heat 
ratio is 7 = 5/3, and the cooling times are t coo \£l — 3 (top 
left), tcooifi = 5 (top right), t coo \Sl = 6 (bottom left), and 
tcooifi = 7 (bottom right). In this particular case we did 
not complete a t coo \£l = 4 run since the icooi^ = 5 run 
had already shown signs of fragmentation prior to the com- 
pletion of the tcooifi = 4 simulation. For tcooi^ = 3 there 
are a large number of fragments, consistent with lRice et alJ 
(2003). For t coo iQ = 5 there are a number of fragments, 
while for t coo i£l = 6 there is only a single fragment that 
in the image can be seen just below the central object. For 
icooifi = 7, which ran almost 7 outer rotation periods, the 
disc is clearly unstable at all radii, but there is no signs of 
fragmentation. 
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—9.7 x 10~ 7 


0.1 


5/3 
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— 1.0 x 10 — 7 
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—3.8 x 10~ 5 
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5/3 
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no clumps 


0.1 


7/5 
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—8.8 x 10~ 7 
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7/5 
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—6.6 x 10~ 8 
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13 


no clumps 


0.25 


5/3 
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-9.4 x 10~ 6 


0.25 


5/3 
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-3.0 x 10~ 7 


0.25 


5/3 


7 


no clumps 


0.25 


7/5 


11 


-8.2 x 10~ 7 


0.25 


7/5 


12 


-7.2 x 10~ 7 


0.25 


7/5 


13 


no clumps 


0.5 


5/3 


6 


-4.9 x 10~ 5 


0.5 


5/3 


7 


no clumps 


0.5 


7/5 


11 


-1.0 X 10~ 5 


0.5 


7/5 


12 


-7.5 x 10~ 6 


0.5 


7/5 


13 


no clumps 



Table 1. Table showing the results of a series of simulations con- 
sidering discs with masses between iWdisc = 0-1 an d A^disc = 0.5, 
specific heat ratio of 7 = 5/3 and 7 = 7/5, and various cool- 
ing times. These results suggest that the fragmentation boundary 
does not depend on disc mass, and that for 7 = 7/5 fragmenta- 
tion may occur for cooling times almost twice the local dynamical 
time. 



We repeated the above procedure for disc masses of 
Mdi S c = 0.25, and Mdi sc = 0.5 and for specific heat ra- 
tios of 7 = 5/3 and 7 = 7/5. The results are summarised 
in Table 1. The columns in Table 1 are the ratio of disc to 
central object mass, Mdisc/M,, the specific heat ratio, 7, the 
cooling time, t coo \Q, and if fragmentation occurs, the total 
energy (in code units) of the densest clump, Eu>t , where Eu>t 
is the s um of the thermal energy and gravitation al potential 
energy llBate et alJll995Jl . In the earlier work of lRice et alJ 
( 2003) there was a suggestion that the cooling time required 
for fragmentation may depend on the total disc mass, rela- 
tive to the mass of the central object. The results shown in 
Table 1 suggest that there is no disc mass dependence. Frag- 
mentation occurs for t C ooi^ between 6 and 7 when 7 = 5/3 
and between 12 and 13 when 7 = 7/5, for all disc masses 
considered . The reason why there is a difference between 
IRice et all ll2003l) is unclear. Their discs had slightly steeper 
surface density profiles (E oc r -7 ^ 4 rather than E oc r" 1 ), 
and it is possible that their t coo \il = 5 simulation that did 
not fragment, may have done so had it been run for longer. 
Table 1 also shows that in all the cases where clumps were 
detected, the total energy of the densest clump is negative 
and that at least the densest clump is bound. 

Although Table 1 shows that the fragmentation bound- 
ary occu r s for cooling times longer than that predicted by 
iGammiel l|2001) , for 7 = 5/3 the required cooling time is 
still smaller than the local dynamical time. It also shows 
that as the specific heat ratio decreases, the required cool- 
ing time increases and is almost twice the local dynamical 
time for 7 = 7/5. The fragmentation boundary for a disc 
mass of Mdisc = 0.25 and for both of the specific heat ratios 
considered is shown in Figures |5] and |3] Figure shows the 
final surface density structures for Mdi sc = 0.25, a specific 
heat ratio of 7 = 5/3, and for cooling times of t coo i£l = 6 
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Figure 2. Surface density structure of discs with a mass of 
A^diac = 0.25, a specific heat ratio of 7 = 5/3, and cooling 
times of t coo in = 6 (left hand panel) and t coo iQ = 7 (right hand 
panel). The lack of fragmentation in the right hand panel suggests 
that the fragmentation boundary is at a cooling time of between 
*cool^ = 6 and i coo iQ = 7. The colour scale of the density and 
the linear scale of the image are the same as in Fig. 1. 




Figure 3. Surface density structure of discs with a the same mass 
as in Figure |3 but with a specific heat ratio of 7 = 7/5, and 
cooling times of t coo \Q = 12 (left hand panel) and t coo iQ = 13 
(right hand panel). The lack of fragmentation in the right hand 
panel in this case suggests that for 7 = 7/5 the fragmentation 
boundary is at a cooling time of between t coo iQ = 12 and t CO ol^ = 
13. The colour scale of the density and the linear scale of the image 
are the same as in Fig. 1. 



(left panel) and t coo i^ = 7 (right panel). The t coo \£l — 7 
simulation was evolved for almost an outer rotation period 
longer than the icooifi = 6 simulation yet shows no signs of 
fragmentation. The discs shown in Figure [3] have the same 
parameters as in Figure [5] except 7 = 7/5, and the cooling 
times are t CO oiH = 12, and t C ooi£2 = 13. Again there is no 
sign of fragmentation in the right hand panel which was also 
evolved for almost an outer rotation period longer than the 
simulation shown in the left hand panel. 

As a further numerical check, we repeated one set of 
calculations using 125000 particles rather than 250000 par- 
ticles. We considered only the case where Mdisc = 0.25 and 
7 = 5/3. The result with 125000 particles was the same as 
the simulation with 250000 particles. Fragmentation occured 
for f C ooi = 6f2~ and did not occur for t coo i = 7fi . There- 
fore, not only do the simulations that fragment satisfy the 
Jeans Criterion for fragmentation feate fc Burkertll997l) . it 
appears that the results are resolution independent. 




Figure 4. The relationship defined by equation^for 7 = 2 (solid 
line), 7 = 5/3 (short-dashed line) and 7 = 7/5 (long-dashed line). 
The data points show the couples (tf rag f2, « max ) as derived from 
the simulations: the green squares refer to our simulations, while 
the blue triangle refers to iGammiel J2001J). The horizontal green 
line illustrates the constant values a = 0.06251. 



3 A MAXIMUM VALUE FOR 
GRAVITATIONAL STRESSES 

Based on the results summarized in Table for every 
value of Mdisc/Mt and 7, we can define a minimum cool- 
ing time for which no fragmentation occurs, i n f and a max- 
imum cooling time for which fragmentation does occur, tf. 
The boundary value of t coo i for fragmentation can there- 
fore be denned as tf rag = l/2(t n f + if), with a correspond- 
ing uncertainty given by Aif rag = l/2(t n f — if). The stress 
corresponding to if rag , can be computed from equa- 
tion and the corresponding uncertainty is given by 
(amaxAf rag )At frag . The resulting values of t frag 
and a max are shown as data points in Fig. [I] together with 
the curves defined by equation Q, for three values of 7 =2, 
5/3 and 7/5. The filled green squares with error bars refer to 
the simulations presented here. The open blue triangle rep- 
resent the value found by IGammiel i200 ll) in his local, 2D 
simulations that assumed 7 = 2. This is consistent with our 
result which suggests that for 7 = 2, fragmentation should 
occur between t coo iQ — 3 and tcooiQ = 4. In fact, it is worth 
noting that since Gammie's simulations are 2D, we should 
not expect a perfect agreement between our 3D results and 
his ones. This can be partially seen already from Fig. ^] In 
particular, care should be taken in considering the role of the 
adiabatic index 7, which has a different physical interpreta- 
ti on in 2D and in 3D. However, as discussed in more detail 
in IGammiel <l200ll) . a mapping is possible between the 2D 
and the 3D adiabatic indeces. In the case of self-gravitating 
discs, Gammie's choice of a 2D adia batic index equ al to 2 
does correspond to 7 = 2 also in 3D dGammidbOOlh . 

As can be seen, fragmentation occurs at an almost con- 
stant value of a (a max ~ 0.06, indicated by the horizontal 



Investigating fragmentation in discs 5 



green line in Fig. , thus vindicating the idea that gravita- 
tional instabilities cannot provide (in a steady state) a stress 
larger than a max . If the dissipation associated with a max is 
not sufficient to balance the cooling rate, then the reaction 
of the disc is to fragment into bound objects. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we elucidate the processes that lead to the frag- 
mentation of a massive disc. Our main result is the determi- 
nation of a maximum value for the stress that can be pro- 
vided by gravitational instabilities in a quasi-steady state. 
We then argue that fragmentation will occur whenever the 
external cooling requires, in order to be balanced by internal 
heating, a stress larger than this maximum value, that we es- 
timate to be Qmax ~ 0.06 (in units of the local disc pressure). 
As a consequence, discs with larger values of the ratio of the 
specific heats will be less susceptible to fragmentation. For 
7 = 7/5, for example, we estimate the fragmentation cooling 
time to be between 11S7 -1 and 12Q~ , com pared to between 
and AO,' 1 for 7 = 2 <Gammidl200lh . 
We wish to stress that the threshold value for a that 
we have found here refers to a quasi-steady state, in which 
the disc stays in thermal equilibrium and the relevant phys- 
ical quantities do not change significantly on time scales 
shorter than the ther mal timescale. We have already shown 
jLodato fc Ricell2005l) how very massive discs (with masses 
comparable to that of the central object) can generate tran- 
sient strong spiral episodes, with correspondingly large val- 
ues of the stress a, which, however, do not last for longer 
than one dynamical timescale (see details in lLodato fc Ricel 
120051) . 

A further remark is in order, in reference to the pos- 
sibility of non-local transport in self-gravitating discs. In 
all our simulations, we did not find any significant evi- 
dence for non-local t r anspo rt of energy due to self-gravity 
jLodato fc Ricdl2004l200Sh . If the disc does not fragment, 
the dissipation provided by the gravitational stresses bal- 
ances almost exactly the imposed cooling rate. However, this 
conclusion might depend on the simulation setup and, in 
parti cular, on the assum ed radial dependence of the cooling 
time. iMei ia et alJ (120051) claim to find evidence for non-local 
energy transport in their simulations that employ a t coo \ con- 
stant with radius (rather that oc fi -1 , as we do). If this is 
the case, it might offer a possible escape route for the disc 
in order to avoid fragmentation. Consider the case where 
icooi^ is such that the disc is stable against fragmentation 
in the inner regions, but would fragment in the outer regions, 
following our prediction in equation (which is based on 
the implicit assumption that energy dissipation is viscous). 
In such a situation, the inner disc could "help" the outer 
disc, by heating it up via non-local effects and preventing 
fragmentation. This process is not viable in the simulations 
presented here, since here the whole disc is uniformly sta- 
ble or unstable with respect to fragmentation and increasing 
the cooling rate of the inner disc via non-local energy trans- 
port would make it fragment. Clearly, further simulations of 
the fragmentation process in discs with a radius-dependent 
icooifi are needed to investigate this issue. 

Finally, w e note that th e fragmentation requirements 
determined bv lRafikovl j2005tl were calculated assuming the 



much shorter cooling times of iGammiel i200lt) . and by as- 
suming that the Toomre Q parameter must be unity or less 
for fragmentation to commence. Since the actual Q value 
required for fragme ntation can be slightly greater than 1 
JPickett et alj ll998). and since t he required coolin g time in 
discs with 7 = 7/5 (as used bv iBosd (Il998t 120021) ') can be 
larger than that predicted by iGammiel 1200 J) , the condi- 
tions for fragmentation may not be as stringent as those 
determined bv lRafikovl (2005). 
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